function mf = QLaw_MEE_fun(input, auxdata, sc, st, WP, r_pmin, del_cond)

%% gravitational constant and spacecraft 
mu      = auxdata.mu;    % km^3/s^2
Re      = auxdata.Re;              % km
Tmax    = auxdata.Tmax;                   % N
Mmax    = auxdata.Mmax;                 % kg
Isp     = auxdata.Isp;                % sec
g0      = auxdata.g0;                % m/s^2
c       = auxdata.c;
F       = auxdata.F;
AU2km   = auxdata.AU2km;

%% Weight definition
% WP = 1;
Wa = input(1);
Wf = input(2);
Wg = input(3);
Wh = input(4);
Wk = input(5);

W = [WP; Wa; Wf; Wg; Wh; Wk];

%% scaling parameter
k = input(6);
m = input(7);
n = input(8);
r = input(9);

params = [k; m; n; r; r_pmin];

dL = 10 * pi/180;

i = 1;

while 1

    F = Tmax / sc(7);
    i = i+1;
    %% Q
    Q = proximity_quotient(sc, st, W, params, F, mu);
    
    %% Q gradient and control definition
    [dQdt, u_opt] = proximity_quotient_dot(sc, st, W, params, F, mu);

    uu = u_opt / norm(u_opt) * F;
    dQdt = dQdt / norm(u_opt) * F;
    
    %% dt definition
    pp = sc(1);
    ff = sc(2);
    gg = sc(3);
    hh = sc(4);
    kk = sc(5);
    LL = sc(6);
    
    [~,ecc,~,~,~,nu] = lowThrustMee2Coe(pp,ff,gg,hh,kk,LL);
    
    rr = pp / ( 1 + ecc * cos(nu) );  % radius from the central body
    dt = dL * sqrt(rr^3/mu);

    if Q + dQdt*dt < 0
        dt = - Q / dQdt;
    end
    
    k1 = x_dot(sc, mu, uu, Tmax, c);
    k2 = x_dot(sc + k1*0.5 * dt, mu, uu, Tmax, c);
    k3 = x_dot(sc + k2*0.5 * dt, mu, uu, Tmax, c);
    k4 = x_dot(sc + k3 * dt, mu, uu, Tmax, c);

    sc = sc + dt*(k1 + 2*(k2 + k3) + k4)/6;
    sc = real(sc);
    
    %% target trajectory
    k1 = x_dot(st, mu, [0;0;0], 0, 0);
    k2 = x_dot(st + k1*0.5 * dt, mu, [0;0;0], 0, 0);
    k3 = x_dot(st + k2*0.5 * dt, mu, [0;0;0], 0, 0);
    k4 = x_dot(st + k3 * dt, mu, [0;0;0], 0, 0);

    st = st + dt*(k1 + 2*(k2 + k3) + k4)/6;
    st = real(st);
    
    %% break condition check    
    cond = abs(sc-st);
    if sum(cond(1:5) < del_cond) == 5
        break
    end
    
    if sc(end) < 0
        sc(end) = 0;
        break
    end
    if i > 100000
        sc(end) = 0;
        break
    end
end

if sc(end) >= 1
    mf = 0
elseif sc(end) >= 0
    mf = -sc(7,end)
else
    mf = 0
end